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Abstract 

KADATH is a library that implements spectral methods in a very modular 
manner. It is designed to solve a wide class of problems that arise in the con- 
text of theoretical physics. Several types of coordinates are implemented and 
additional geometries can be easily encoded. Partial differential equations of 
various types are discretized by means of spectral methods. The resulting 
system is solved using a Newton-Raphson iteration. Doing so, KADATH is able 
to deal with strongly non-linear situations. The algorithms are validated by 
applying the library to four different problems of contemporary physics, in 
the fields of gauge field theory and general relativity. 
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1. Introduction 

Spectral methods are a class of numerical methods that aim at solving 
partial differential equations. For detailed presentations of those techniques, 
the reader should refer to the numerous books available like [i], [sl, 0, HI, [oj. 
The basic idea is to describe any field by an appropriate linear combination of 
known functions called the basis functions. Classic examples of basis are the 
trigonometrical functions or orthogonal polynomials (Chebyshev, Legendre, 
etc.). The description of functions by their spectral expansion is by essence 
non-local, which is to be contrasted with finite difference schemes. Spectral 
methods are particularly appealing because of the very fast convergence of 
the spectral expansion to the real function it describes. More precisely for 
C°°functions the error decreases faster than any power-law of A^, where N 
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is the order of the expansion. A mult i- domain decomposition, where the 
physical space is decomposed into several computational regions, is usually 
used to ensure such smoothness of the functions. 

Originally, spectral methods were used in the context of numerical hy- 
drodynamics and led to the successful computations of turbulence regimes 
and various two or three-dimensional flows. Numerous applications can be 
found in [2I, 0]. The application of spectral methods to the field of general 
relativity is somewhat more recent and starts with the pioneer work of the 
Meudon group, in the late eighties. Since then, such methods have been 
successfully applied by several groups to systems like binary black holes or 
magnetized neutron stars. A review of spectral methods in the context of 
general relativity can be found in Q. 

LORENE is the library written by the Meudon group that implements spec- 
tral methods for general relativity. It is written in C++, is used by several 
groups worldwide and produced many results (see [sl for references). Nev- 
ertheless, with time, it appeared that LORENE did encounter some difficulties 
and had some serious limitations. The first one comes from the fact that 
the library is designed to work with spherical coordinates only and from the 
difficulty to implement new geometries like the bispherical one for instance. 
Second, LORENE solves systems by an iterative loop that requires a splitting 
of the equations in terms of operators (like the Laplace one) and sources. 
This splitting is obviously not unique. If it can be natural in some cases, this 
is not true when strong non-linearities occur, like for gauge field problems 
where the use of LORENE proved problematic (the example shown in Sec. 17.11 
was never successfully computed by LORENE). Finally, and even if learning 
LORENE is manageable, it can still be somewhat difficult to write a complete 
code. A more user-friendly coding style would be a good thing. 

For all those reasons, it has been decided to think about a new library 
that would use the many successes of LORENE and try to improve on its 
weaknesses. The KADATH library is the result of this process and this work 
is the first paper devoted to its description. The design of KADATH was also 
inspired by some aspects of other spectral solvers, like the ones described in 



lOl . Three different types of coordinates have been coded so far but 
the setting is such that the inclusion of additional geometries is relatively 
easy and straightforward. Spectral expansion is done either with respect to 
Chebyshev or Legendre polynomials and trigonometrical functions. Systems 
of equations are solved as such, by means of a Newton- Raphson iteration and 
do require only a minimal rewriting of the equations. They are passed to the 
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solver as character strings. Those strings are interpreted by KADATH using a 
syntax close to the LateX one [lij] , so that the interface should be relatively 
clear to any modern physicist. The library is written in C++. It has been 
made public and can be downloaded by going to the KADATH web-site 13 . 



In its current state of development, the library deals with boundary value 
problems. Moreover, the geometry of the various boundaries must be known 
in advance and is fixed during the whole computation. Such kind of setting is 
useful in computing stationary or periodic solutions and encompasses a wide 
class of problems, as the four examples discussed in Sec. [7] illustrate. Never- 
theless, the field of application of KADATH would be significantly broadened if 
some of those restrictions could be lifted. One possible extensions would be 
the possibility to deal with free boundary problems where the geometry is no 
longer fixed but numerically determined in the course of the computation. 
This would be required, for instance, to compute neutron stars in rotation, 
where the shape of the surface of the star is an unknown. Another addition 
to KADATH would be the inclusion of tools to study time-evolution problems. 
This is however a major task that would require much time and work. Those 
possible extensions, as long as a few others, are discussed in more detail in 
Sec. El 

This paper contains two main parts. In the first one, the basic numerical 
techniques used by the library are exposed. Not all the details are given, for it 
would be cumbersome. Nevertheless it should give the reader a good feeling 
of what KADATH is about. Section [2] describes the way multi-domain settings 
are implemented along with the three different types of geometries currently 
present in KADATH . In Sec. [3l after a short introduction on spectral expansion, 
the various basis used by KADATH are detailed, for different geometries and 
different types of mathematical objects (scalars, tensors, etc.). The way 
additional regularities are enforced is also discussed in some details. Section 
m is devoted to the discretization of partial differential equations by r and 
Galerkin methods. The resulting non-linear system is solved by a Newton- 
Raphson iteration which is described in Sec. [51 The computation of the 
Jacobian and its inversion are discussed, especially in the context of parallel 
computing. Section El gives an outline of what a code that uses KADATH should 
look like. 

The second part of the paper is concerned with four different test problems 
(Sec. [7]). Those problems have been chosen because they illustrate different 
aspects of KADATH , in terms of the type of equations, variables or geometries. 
If those are not exactly new results, they are far from being trivial and 
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they all relate to very contemporary physics. The existence of vortons in a 
particular gauge field theory is confirmed in I7.H a critical collapse solution 
is constructed in Sec. \7.2\ binary black holes spacetimes are obtained in 17.31 
and the Kerr solution for a single rotating black hole is computed in 17.41 
Future developments and applications are discussed in Sec. [81 

2. Setting the geometry 

2.1. Multi-domain decomposition 

KADATH implements spectral methods in a multi-domain manner where the 
physical space is split into several computational domains. The advantages 
of using a multi-domain setting are numerous. First of all, when dealing 
with discontinuous fields, the domains can be chosen so that the surfaces of 
discontinuities coincide with the boundaries. Doing so, in each domains, all 
the functions are C°°, thus recovering the spectral convergence that would 
be lost otherwise (see Sec. 13. ip . The use of several domains also enables 
to use different resolutions in different parts of space, increasing accuracy in 
regions where needed. This is called fixed- mesh refinement. There are also 
some cases where setting a global and regular set of numerical coordinates is 
troublesome. This is for instance the case of the bispherical coordinates, as 
can be seen in Sec. 12.31 Such difficulties can be overcome by setting different 
numerical coordinates in different regions of space. 

For each domain, there is a mapping between a set of physical coordinates 
(the (r, 9,ip) of spherical coordinates for instance) and a set of numerical co- 
ordinates. The spectral expansion is performed with respect to the numerical 
ones. The mapping between the two sets of coordinates ensures that the nu- 
merical ones lies in the appropriate range to do the expansion. For instance, 
a coordinate must be in [—1, 1] if Chebyshev polynomials are used. 

The physical space is divided into a finite set of domains. Only touching 
and not overlapping domains are considered. The domains are described by 
the class Domain and its derived classes. The way the various domains are 
set with respect to each other is encoded in the abstract class Space and 
its derived classes. This is needed, for instance, to write the appropriate 
matching conditions across the boundaries of the domains. 

KADATH provides several domain decompositions that are listed below. 
The library has been designed to be very modular in terms of the geometry, 
so that it is relatively easy to implement new types of domains and spaces. 
Among the possible future developments, one could think about square-like 
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domains, three-dimensional cylindrical coordinates or spaces with more than 
three dimensions. 



2.2. Spherical space 

The class Space_spheric implements a decomposition of space in terms of 
spherical shells. This is obviously intended mainly for spherical-like objects 
but has also been successfully used for more complicated shapes, like the 
toroidal one (see Sec. I7.ip . This geometry is very similar to the one described 
in H. 

In this setting, the physical coordinates are the standard spherical ones 
{r,6,(f ). One notes 6 the zenith angle and (p the azimuthal one, so that, in 
terms of Cartesian coordinates, one gets: 



For all the type of domains, the numerical angular coordinates {6*, (p*) are 
identical to their physical counterparts. So far, the various fields are supposed 
to be either symmetric or antisymmetric with respect to the z = plane. 
This implies that 6 can be restricted to [0,7r/2]. No symmetry is assumed 
with respect to if, which lies in [0,27r[. Concerning the radial coordinate r, 
three different types of domains are considered. 

• The class Domain_nucleus represents a spherical domain that encom- 
passes the origin of the coordinate system and that extends up to a fi- 
nite radius. The numerical radial coordinate r* relates to r by r = ar*, 
where a is a constant that gives the radius of the domain, r* lies in 



• The class Domain_shell represents a spherical shell, where the radius 
r lies between to finite values. The numerical coordinate r* relates to 
r by an affine law r = ar* + /3, where a and /3 are constants, r* spans 
the intervale [—1, 1]. 

• The last type of spherical domain is called Domain_compact an extends 
from some finite radius up to infinity. This is done by demanding that 



X = r sin 6 cos (p 
y = r sin 6 sin (f 
z = r cos 6 



(1) 
(2) 
(3) 



[0,1]. 



r* relates to r by r 



1 



where a is a negative constant, r* lies 



a (r* — 1) 
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in [—1, 1]. As intended, the domain goes to r = oo, which corresponds 
to r* = 1. 

For axisymmetric problems that do not depend on a variant of the 
spherical space has been implemented. It is called Space_polar and the 
only difference is that the variable if has been removed, making the space 
two-dimensional, thus reducing the computational cost. 

2.3. Bispherical space 

Bispherical coordinates are implemented by the class Space_bispheric. 
In the standard point of view, bispherical coordinates (x, ^, relate to the 
Cartesian ones by: 

asinh?7 , . 

^ = Z (4 

cosh T] — COS X 

a sin X cos ip 



cosh T] — cos X 
a sin X sin ip 



cosh rj — cos X 

The coordinate rj can take all the values in M, whereas x goes from to 
IT. if is the angle around the x-axis and it lies in [0,7r], once the symmetry 
with respect to the plane z = is taken into account. 

Bispherical coordinates are such that the surfaces of constant r] are non- 
concentric spheres located on the x-axis, hence their name. There are given 
by 



^2 

{x — acothrj) + + = t^— (5) 

sinh rj 

This property makes those coordinates very appropriate to deal with physical 
systems consisting of two spherical-like objects. So far, KADATH enables the 
user to consider the space exterior to two spheres (for an application see Sec. 
17.31) . Let us call ri and r2 the radii of those spheres and d their separation. 
It is then easy to see that the space exterior to those two spheres is described 
by the points such that rj G [r/", The values of r/^ and the scale factor a 
are then uniquely defined and obey the following set of equations: 

_ a 

smhr/ = — (6) 
ri 

sinh 77"'" = — (7) 
d = a (cothi]"*" — coth?7~) . (8) 
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Those equations simply translate in terms of rj the positions and sizes of the 
spheres by making use of Eq. Let us note that, doing so, the origin of 
the Cartesian coordinates (i.e. the point x = y = z = 0) cannot be placed 
arbitrarily. 

In bispherical coordinates, spatial infinity is described by r/ = and x = 
implying that the surface r = oo is described by a single line. This makes 
a simple compactification like the one used in the spherical case (see Sec. 
12.21) impracticable. This can be seen if one forms the ratio of the Cartesian 
coordinates with respect to r. As one approaches spatial infinity such ratio 
goes like 

-^^^ (9) 

which does not admit a well-defined limit when rj ^ and x ~^ 0. This 
effect is closely related to the fact that the coordinate transformation (jlj) si 



only C° and not C°°at infinity, as can be seen in H 



One way to overcome this difficulty is to excise a region of space around 



the point rj = and X = (see ll|] for another possibility). Doing so. 



the physical space does no longer extend up to infinity but only to a finite 
distance. The shape of the newly created exterior boundary is defined by the 
shape of the excision region. For instance, in order to get a spherical outer 
boundary of radius R, one needs to excise the region: 



sinh^ 1] + sin^ x 



2 ^ (10) 

(cosh r] — cos x) ^ 

The presence of the exterior boundary can also be useful for some problems 
for which outer boundary conditions must be prescribed. This is especially 
true for problems with radiation where outgoing waves can be present. 

To summarize, KADATH implements bispherical coordinates in the region 
defined by r^^ < 77 < 77+ and Eq. ( ITOj) . This corresponds to the region of 
space exterior to two spheres of radii r\ and r2 and distant from d and inside 
an outer sphere of radius R. The center of the outer sphere coincides with the 
origin of the coordinates and is not freely specifiable. However, if needed, Eq. 
( fTOj) could be modified to accommodate any smooth outer boundary. Both 
the inside of the two inner spheres and the exterior of the outer one, can be 
matched with spherical-like domains as those seen in Sec. l2.2l to describe the 
desired physical space. 
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In order to perform numerical expansions, one needs to map the bispher- 
ical coordinates {t],x) onto some numerical coordinates {ri*,x*)- Unfortu- 
nately, with the region previously defined, this cannot be done simply using 
a single domain. The bispherical region must be split into several domains, 
each of them being mapped onto squares of numerical coordinates. In order 
to do the splitting the following types of domains are defined. 

• The class Domain_bispheric_rect represents a rectangular domain in 
terms of [rj, x)- More precisely, 1] lies in between two finite values r^mm 
and ?7max- It relates to the numerical coordinate rj* by an affine-law 

V = ^ V + ^ (11) 

so that r]* G [—1, 1]. x S^es from a lower boundary Xmin up to n and 
relates to its numerical counterpart via 

X = (Xmin - vr) X* + vr (12) 

so that X* ^ [0)1]- This choice of mapping will be justified by the 
chosen basis of decomposition in Sec. I3.2.2I 

• The class Domain_bispheric_chi_f irst describes domains for which 
T] is given as a function of x- X S^es from to a boundary Xmax and 
relates to x* by 

X = XmaxX* (13) 

X* lies in [0, 1]. r] goes from fixed value r/iim up to a variable bound 
/ (x). The function / (x) is given by the surface of excision ( ITOl) and 
is computed numerically, t]* goes from [—1,1] and is given by 

^ ^ /(x)-^lim ^. /(x)+^lim ^^^^ 

• In the class Domain_bispheric_eta_f irst, the variable x is given in 
terms of rj. r] is located between ?7min and r^max and relates to rj* simply 
by Eq. ( ITTj) . x S^es a variable bound g {rj) up to n. As for the 
domain Domain_bispheric_chi_f irst, the function g (77) is defined by 
the excision shape. One then gets 

X = {9 iv) - vr) x'' + TT (15) 

with X* e [0, 1] 
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Table 1: Relation between the various domains of the bispherical space. The table shows 
to what correspond the boundaries of each domain. The labels are the same as in Fig. [T] 



Domain 


rj* = —1 


r]* = 1 


x* = o 


x' = i 


1 


left inner sphere 


outer sphere 


X-axis 


domain 2 


2 


left inner sphere 


domain 3 


X-axis 


domain 1 


3 


domain 2 


domain 4 


X-axis 


outer sphere 


4 


right inner sphere 


domain 3 


X-axis 


domain 5 


5 


right inner sphere 


outer sphere 


X-axis 


domain 4 



In the three types of domains (p coincides with the associated numerical 
coordinates and goes from to vr. 

A set of five such domains is needed to describe the space between the 
two inner spheres and the exterior one. This is done in the following way. 
Let us pick two sets of values (?7i,Xi) and (^725X2) that lie on the excision 
region. This means that they correspond to circles on the exterior sphere. 
One will assume that r/i < and that 772 > 0. The domain decomposition is 
as follows, the labels of the domains are given in Fig. [1) 

• One Domain_bispheric_chi_f irst for which < % < Xi (domain 1). 

• One Domain_bispheric_rect for which Xi ^ X ^ ^ and f]~ < f] < f]i 
(domain 2). 

• One Domain_bispheric_eta_f irst for which r]i < t] < r]2 (domain 3). 

• One Domain_bispheric_rect for which X2 < X ^ ^ and < r] < 7]^ 
(domain 4). 

• One Domain_bispheric_chi_f irst for which < x ^ X2 (domain 5). 

As already stated, spherical domains can be added to this decomposition to 
account for the interior of the inner spheres or for the exterior of the outer 
one. This setting is shown in Fig. [T]both in the (r/, x) plane and in terms of 
Cartesian coordinates. Table [1] shows to what correspond the boundaries of 
each domain. In particular, one can see that the mappings have be chosen 
so that X* = always corresponds to the x-axis. This property will prove 
useful when regularity conditions will be enforced (see Sec. I3.4.3p . 
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Figure 1: Multi-domain decomposition of the bispherical space. The inner spheres have 
radii of 0.5 and 1 and are separated by a distance of 2. The radius of the outer sphere is 
set to 4. The left panel shows the bispherical coordinates, along with the excision region 
around 77 = and x = 0. The associated Cartesian coordinates, in the {x, y) plane, are 
shown on the right panel. 

2.4- Cylindrical space 

Standard three-dimensional cylindrical coordinates are currently not im- 
plemented in KADATH . The setting exposed in this section is devoted to the 
study of critical phenomena in general relativity (see Sec. 17.21 and references 
therein). The sought critical solution can be described on a two-dimensional 
space of cylindrical topology. If this setting is more specialized than spherical 
or bispherical coordinates, it is presented here as an illustration of the ability 
of KADATH to deal with various geometries. 

The cylindrical space is described by the class Space_critic. The two 
coordinates are x that ranges from to 1 and r that goes from to 27r. x 
relates to the height of the cylinder and r is the azimuthal angle. Near x = 
the fields of interest are either symmetric or antisymmetric. In order to take 
this symmetry into account, space is split into two types of domains: 

• Domain_critic_inner where x lies in [0,Xiim] and simply relates to the 
numerical coordinate x* by x = xumX*. x* covers the interval [0, 1]. 

• Domain_critic_outer where x is in [xum, !]• The numerical coordinate 
X relates to x by an arfane law x = x H and spans 
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[-1,1]- 



In both domains, the numerical angular coordinate r* coincides with r. 
Given the fact that the fields have some symmetries with respect to r = tt, 
it is possible to restrict r* to [0, 7r[. 

3. Spectral expansion 
3.1. Generalities 

An extensive discussion of spectral methods is beyond the scope of this 
work. Only the basic properties required for the comprehension of the paper 
are presented. Refs. [H, y, 0, 0, ^ ^ could be useful for the reader interested 
in more details. 

In the one-dimensional case, let / (x) be a function on an interval A. 
Given a set of known orthogonal functions $j (x) on A, spectral theory en- 
ables to construct an approximation of / in terms of a finite sum of the 
$j (x). This sum is called the interpolant of / and is expressed as: 

N 

lNf{x) = Y,^^^d^) (16) 

where is the order of the approximation. Typically, the basis functions are 
either orthogonal polynomials like Legendre or Chebyshev ones, or trigono- 
metrical functions. In this latter case, spectral theory is nothing but the 
theory of discrete Fourier transform. 

It can be shown that there exist A^ + 1 points Xi inside A such that / and 
its interpolant exactly coincide at those points 

/(x,)=J^/(x,). (17) 

Those points are called the collocation points. In KADATH , one works with 
the so-called Gauss-Lobato points which ensure that the boundaries of the 
interval A coincide with the first Xq and last x^v collocation points. 

Spectral theory provides a rule to compute the coefficients a, in terms 
of the values at the collocation points. Thanks to Eq. (1161) . the reverse 
operation, that is the computation of the / (xj) in terms of the a^, is also 
possible. So a function / can be described either in terms of its coefficients 
or in terms of its value at the collocation points. Depending on the operation 
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to be performed on /, a description could be more useful that the other. One 
is working in the coefficient space when using the Oj and in the configura- 
tion space when considering the f {xi). The two descriptions are completely 
equivalent in the sense that there is no information lost when going from one 
space to the other. 

The most appealing feature of spectral methods is the very fast conver- 
gence of the interpolant I^f to the real function /. Indeed, it can be shown 
that, if / is a C°°function, then I^f converges to / faster than any power-law 
of N. This is known as the spectral convergence. In the case of analytic func- 
tions, the convergence is even exponential. Such fast convergence enables to 
reach good accuracy with only moderate number of points, especially com- 
pared to other methods like finite difference schemes. As already stated in 
Sec. 12.11 an appropriate multi-domain setting can usually ensure that the 
functions are C°°in each domain even for globally less regular functions. Fail- 
ing to do so will make the convergence only follow a power-law. This is called 
the Gibbs phenomenom in the case of a discontinuous function. 

3.2. Scalar fields 

The choice of basis functions is crucial to the success of any spectral 
solver. Moreover it is often a very good way of checking equations. Indeed, 
for complicated ones, where many terms are involved, all of them must end 
up having consistent basis. This section is devoted to the case of scalar fields. 
Higher order tensors are discussed in Sec. 13. 3[ 

3.2.1. Spherical coordinates 

In the case of spherical coordinates, the standard basis of decomposition 
of a scalar field is obtained by demanding that this field can be expressed as 



a polynomial in terms of the Cartesian coordinates, as what is done in [15 
This condition prevents the appearance of singularities on the axis and at 
the origin. In addition, one requires that the fields are either symmetric or 
antisymmetric with respect to the plane z = 0. This last assumption covers 
most of the situations of interest. By expressing the Cartesian polynomials 
in terms of the spherical coordinates, one can get some constraints on the 
basis of expansion. 

Given that no symmetry is assumed with respect to the azimuthal angle 
(f, a standard discrete Fourier transform is performed. The basis of decompo- 
sition consists of the trigonometrical functions cos (rrup) and sin (rrvp) . The 

collocation points are the o9j = 27i with < i < A^. 
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Concerning the zenith angle 6, the symmetry with respect to the plane 

2; = is taken into account. This implies that the collocation points are the 
vr j 

6j = -^j^- The associated basis functions are trigonometrical functions of 

one given type and of one given parity. The exact choice depends on both 
the function and the (p basis. More precisely, for a symmetric function, one 
uses even cosines cos (2j^) when the basis in ip is such that m is even and 
odd sines sin ((2j + 1)6') otherwise. For antisymmetric functions, odd cosines 
are used for m even and even sines otherwise. 

The basis with respect to the radial numerical coordinate r* depends on 
the type of domain. For both the shells and the compactified domain, it is 
the standard Legendre of Chebyshev polynomials with the associated Gauss- 
Lobato collocation points. In the nucleus, only polynomials of a given parity 
are considered, in order to impose some regularity at the origin. If the asso- 
ciated basis with respect to 6 is even (resp. odd), with either sines or cosines, 
even (resp. odd) polynomials are used. This is true for both symmetric and 
antisymmetric functions. The collocation points are the Gauss-Lobato points 
that are located in [0, 1]. 

Let us note that with this choice of basis, not any function is a true 
polynomial of the Cartesian coordinates. The choice made here is slightly 
less restrictive and has the advantage of being convenient and easy to handle. 
Additional regularity conditions are discussed in Sec. 13.41 

3.2.2. Bispherical coordinates 

The same guidelines as in the spherical case are used to derive appropriate 
basis of decomposition. One demands that scalar fields can be expressed as 
polynomials in terms of Cartesian coordinates. Once again, only symmetric 
or antisymmetric functions with respect to the plane z = are considered. 
This gives the following constraints on the basis. 

The angle if is expanded on cosines cos {rrvp) (resp. sines) for symmetric 
(resp. antisymmetric) functions. The associated collocation points are the 

= 4. 

For the coordinate x*^ due to regularity conditions on the x-axis, only 
polynomials (Chebyshev or Legendre) of some given parity must be taken 
into account. The parity depends on the basis with respect to ip. Even 
polynomials must be used when m is even and odd otherwise. This can be 
seen in the expressions of the Cartesian coordinates given by Eqs. (jlj). This 
is true for both symmetric and antisymmetric functions. The collocation 
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points are the Gauss-Lobato points that are located in [0, 1]. 

The coordinate t]* is expanded on standard polynomials, either Chebyshev 
or Legendre. The collocation points are the usual Gauss-Lobato points and 
span the interval [—1, 1]. 

As in the spherical case, this choice of basis is slightly less restrictive than 
what would be needed to get true Cartesian polynomials (see Sec. 13.41 for 
more regularity conditions). 

3.2.3. Cylindrical coordinates 

As already stated, the functions considered in the critical phenomenom 
case (see Sec. I7.2p are either symmetric or antisymmetric with respect to 
X = 0. To account for this fact, in the inner domain, only polynomials of the 
appropriate parity are used, either Chebyshev or Legendre. The collocation 
points in terms of x* are the Gauss-Lobato points and are located in [0, 1]. 
In the outer domain, standard polynomials are used. 

Some symmetry is also taken into account with respect to r = vr. If 
a function / is such that / (r + vr) = / (r) then only even trigonometrical 
functions are considered and / is expanded on both cos (2zr) and sin(2zr). 
The other possibility is a function g such that g {t + tc) = —g (r). In that 
case, only odd trigonometrical functions are used. 

3.3. Higher rank tensors 

For tensors, the choice of spectral basis depends on both the geometry 
and the tensorial basis. The simplest choice is the one of a Cartesian tensorial 
basis. Let us precise that this choice can be made almost independently of the 
geometry, as long as Cartesian coordinates are defined. For instance, a vector 

— * 

V can be described by its Cartesian components {V^ ,V^) in spherical 
geometry, meaning that each component is given in terms of (r, 6, (p). 

Let us first turn to the case of a vector with a Cartesian tensorial basis. 
The appropriate basis are obtained by demanding that this vector can be 
expressed as the gradient of a scalar field. Given this assumption, it simply 
follows the each component can be expressed as a polynomial of the Cartesian 
coordinates. Thus each component behaves like a scalar field and the same 
spectral basis are used. As far as the symmetry 2; = is concerned (for the 
spherical and bispherical cases), it is assumed that the components x and 
y of the vector are symmetric and the z one antisjTiimetric (i.e. the vector 
is the gradient of a symmetric scalar field). The same is true for higher 
order tensors in Cartesian tensorial coordinates for they can be obtained as 
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tensorial product of vectors. When needed, the appropriate z = symmetry 
is also obtained from this fact. 

The case of a spherical orthonormal tensorial basis is also implemented, 
when one is working in spherical geometry. The appropriate spectral basis of 
each component can be derived from the Cartesian case by making a careful 
use of the passage formulae that relates the two basis. For instance, if the 
radial component of a vector behaves like a (symmetric) scalar, this is not 
the case of the 6 one for which the spectral basis with respect to 6 involves 
even sines for m even and odd cosines otherwise, as can be seen from the 
formula 

= V cos ecosip + yy cos esimp- V sin 6. (18) 

In the current state of KADATH , only Cartesian tensorial basis are imple- 
mented in the bispherical case. In the spherical case, both Cartesian and 
orthonormal spherical tensorial basis are defined. The cylindric space only 
deals with scalars so far. Additional cases will be implemented when needed. 

3.4- Additional regularity conditions 
3.4- 1- Galerkin basis 

As stated in Sec. 13. 2[ the spectral basis chosen do not ensure a complete 
regularity of the fields. Some additional constraints must be enforced by 
means of an appropriate Galerkin technique which goes as follows. The fields 
are not expanded onto any set of basis functions but only onto a subspace, 
which verifies the additional constraints one wishes to enforce. For instance, 
let us consider a one-dimensional function / (x) expanded onto even Cheby- 
shev polynomials (x). To fulfill the additional constraint that / vanishes 
at X = it is possible to use the Galerkin basis of the Gi = T2j+2-l- (— 1)*'''^ / 
is expanded onto the Gi and thus, by construction, does verify the constraint 
/ (0) = 0. Let us note that, in general, a Galerkin basis is not orthogonal. 

3.4-2. Regularity on the context of spectral methods 

Given the geometries present in KADATH , two types of regularity must be 
discussed. The regularity on one axis that must be enforced in both spherical 
and bispherical cases and the regularity at the origin r = in the case of 
spherical coordinates only. The main reason why those regularities must be 
carefully handled can be found in the way spectral methods compute some 
ratios. 

As an illustration, let us concentrate on the axis case in spherical coordi- 
nates. For many operators, some functions must be divided by sin 6. This is 



15 



1 r 1 T 1 • r 1 r 1 1-1 C0S6' df 

the case oi the Laplacian oi a scalar / where terms hke -— do appear. 

sm oO 

The division by sin 6* can be troublesome on the axis where it vanishes. In 
the case of spectral methods this difficulty can be overcome by working in 
the coefficient space thus using the fact that the spectral approximation is 
not local. However, doing so, what is computed is not the exact ratio of a 
function g by sin^ (that would diverge for arbitrary function) but the reg- 
ularized one R = - — - — . R is obviously regular for any function g 

smO 

because the finite part of g on the axis has been taken out. It does coincide 
with the real ratio if and only if g vanishes on the axis. 

There is another way to look at this problem. When one computes the 
Laplacian of a function in spectral method, what is actually computed is not 
the real Laplacian, but another operator that includes the finite parts on 
the axis. Those finite parts cause the appearance of additional homogeneous 
solutions (typically solutions that do not vanish on the axis) . Some examples 



of that can be found in [16| , for the radial coordinate, either at the origin or at 
infinity. The extra homogeneous solutions must be dealt with by enforcing 
additional conditions on the solution. If one fails to do so, the resulting 
system will not be invertible. The additional conditions are what we refer as 
regularity conditions and are discussed in more detail in Sec. I3.4.3l and 13.4.41 

3.4-3. Regularity on the axis 

Let us ffist consider the case of a scalar field / in spherical coordinates. 
The regularity requires that / vanishes on the axis, except for the m = 
case. Given the basis of decomposition discussed in Sec. 13.2.11 this can be 
enforced by using the Galerkin basis cos (£6) — 1 instead of standard cosines, 
when m 7^ 0. When the basis with respect to 6 involves sines only, there is 
no need to impose any additional condition. In the case m = the standard 
basis are used. 

For higher order tensors the same guidelines as in Sec. 13.2.11 are used. 
When a Cartesian tensorial basis of decomposition is employed, the regularity 
conditions for each component are the same as in the scalar case. With a 
spherical tensorial basis, they are slightly different. Typically one can allow 
some components (the angular ones) to take non-zero values on the axis 
for higher m. For instance, the component of a vector V relates to the 
Cartesian components by the passage formula f lTSj) . Given the regularity 
conditions for the Cartesian components, it is easy to see that can be 
non-zero on the axis for m = 1. 
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In the bispherical case, the situation is very similar. As seen in Sec. 
12.31 the axis is described by x* = 0. For scalars or Cartesian components 
of tensors, one demands that they vanish on the axis, for m > 0. When 
even Chebyshev (resp. Legendre) polynomials are used, this is done by using 
the Galerkin basis: T2i (x*) + (-1)'^^ (resp. (x*) - (0)). When odd 
functions are used, not additional condition is enforced. 

3.4-4- Regularity at the origin 

The situation at the origin in the case of a spherical geometry is more 
complicated. It can be shown that a true polynomial of Cartesian coordinates 
would vanish as at the origin, where the 6 basis is cos {£6) or sin (£6) (see 
(isl). However this condition is difficult to enforce exactly and, in some cases, 
has been found to generate instabilities. 

As a first step one demands that for £ ^ the functions vanish at the 
origin. When using even Chebyshev (resp. Legendre) polynomials in r, this is 
done by using the Galerkin basis: T2i+2 (r*) + (— 1)*'*'^ (resp. {r*)—L2i (0)). 
Nothing needs to be done if odd polynomials are used. This choice of basis 
enforces the first order of the regularity at the origin. 

However, this is not sufficient to ensure that usual second order operators 
like the Laplacian are well inverted. This is related to the appearance of 
additional homogeneous solutions discussed in Sec. 13.4.21 The second order 
of the regularity condition must be supplemented. This is done by demanding 
that the derivative with respect to r vanishes at the origin, for £ > 1. If this 
is automatically verified for even polynomials, this is not the case for odd 
ones. When dealing with odd Chebyshev polynomials a possible choice of 
Galerkin basis is given by T2j+3 — T2i_^_^ (0)^i) for £ > 1 (and similarly for 
Legendre). 

This choice has proven to be sufficient to ensure that second order differ- 
ential equation are solved properly at the origin. For higher order equations, 
there may be need to further strengthen regularity. Let us note that, con- 
trary to the regularity on the axis, the regularity conditions at the origin are 
the same for scalars or components of tensors. This comes from the fact that 
the transformations that relate the Cartesian components to the spherical 
ones do not involve r (see for instance ( llSp ). 

The regularity on the axis must obviously also be enforced when the origin 
is present. This is resulting in a two dimensional Galerkin basis for r and 6. 
Putting all the pieces together, the appropriate basis of decomposition for a 
scalar field symmetric with respect to z = 0, in a spherical nucleus would be 
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• For Lf : the standard Fourier decomposition in cos {rmp) and sin {rmp) 

• For 9 (regularity on the axis): 

- cos {2je) for m = (£ = 2j) 

- cos {'ijO) — 1 for m even and m 7^ (£ = 2j) 

- sin ((2j + 1) 6) for m odd = 2j + 1) 

• For r (regularity at the origin): 

- T2i (r*) for £ = 

- T^i+i (r^) for £ = 1 

- T2i+2 (r^) + (-1)*"^^ for I even and £ 

- T2i+3 (r^) - (0) Ti for £ odd and £ ^ 1. 

Let us mention that if the regularity on the axis is routinely handled in 
the context of spectral methods, this is not the case of the conditions at 
the origin. For many applications, like the meteorology simulations or the 
black holes with excision (see Sec. 17.31 and 17^ . the origin is not part of the 
computational domain. Even when the origin is present, various tricks are 
usually employed to avoid dealing with regularity conditions. For instance, in 



171 . Il8l ] the region close to the origin is not described by spherical coordinates 
but by Cartesian ones. In [19j, the solution is sought in the configuration 
space with collocation points that avoid the origin. 

3. 5. Effective number of unknowns 

KADATH is designed to solve systems of equations in the coefficient space. 
Using such a setting the unknowns are the true coefficients of the fields of 
interest. By true coefficients one means that some of them are irrelevant. For 
instance when using a standard Fourier transform the coefficient of sin (0) is 
obviously meaningless. Some cases are less trivial. Indeed, the last coefficient 
of an expansion in terms of odd Chebyshev polynomials is also irrelevant. 
This is related to the fact that the value of such a function is, by construction, 
at X = 0. The value of the function is freely specifiable only at — 1 of 
the collocation points. To maintain the bijection between the coefficient 
space and the configuration space one then needs to reduce the number of 
true coefficients. 
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Moreover the regularity conditions discussed in Sec. 13.41 reduce further- 
more the true number of unknowns. As aheady stated, such additional con- 
ditions are imposed by a Galerkin method. This implies that the unknowns 
are not the coefficients onto the standard basis of expansion but rather onto 
the Galerkin basis. The number of coefficients is then different, reflecting 
the fact that the Galerkin basis contains informations about additional con- 
ditions. For instance, suppose one works with + 1 coefficients in terms of 
even Chebyshev polynomials so that the function is expanded onto the 
with < i < A^. Suppose now that the Galerkin basis T2i+2 (r*) + (—1)*'''^ is 
used to enforce that the function vanishes at the origin (like in Sec. I3.4.4p . 
It is easy to see that, in order to keep the same degree of approximation (i.e. 
to truncate the series to the same order in terms of the polynomials), one 
needs to go only up to i = A^ — 1. This is a general feature of the Galerkin 
method. 

The true number of coefficients and so the number of true unknowns is 
determined by KADATH by making use of the type of field considered (scalar, 
vector, rank-2 symmetric tensor...), the symmetry, the tensorial basis and 
the geometry. Such computation is automated and should be transparent to 
the user. 

4. Setting equations 

4.I. The weighted residual method 

The usual way of solving partial differential equations in the context of 
spectral methods is based on the framework of the weighted residual methods 
(see for instance Q for a more detailed presentation). Let us consider an 
equation written formally as i? = 0. R is a general function on the space of 
interest. For instance, if one needs to solve a simple Poisson equation, one 
would have R = Af — S, where 5* is the source and / the unknown field. 
Obviously, more complicated, non-linear cases involving several unknown 
fields are possible. The weighted residual method translates the functional 
equation R = into a finite set of discrete equations by demanding that the 
scalar product (-R, of with respect to some test function ^ vanishes. The 
scalar product is the same as the one used to define the spectral expansion. 

Depending on the choice of test functions one generates different methods. 
KADATH mainly implements the variant known as the r-method. In this case 
the test functions are the same as the spectral basis ones. Suppose that R is 
expanded onto the C, hj R = J2 C C,, C being the coefficient corresponding 
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to the basis function ^. The sum is taken on all the dimensions, up to 
the desired order of approximation. One can then show that the residual 
equations (-R, = is equivalent to demanding that C = 0. Doing so, 
one then obtains as many equations as the total number of basis functions. In 
the r-method, the equations corresponding the last coefficients can be relaxed 
and replaced by equations that describe appropriate matching between the 
domains and boundary conditions. The exact number of equations that must 
be relaxed depends on both the geometry and the equations themselves. 

As mentioned is Sec. 13.41 when regularity conditions must be enforced, a 
variant of the r-method known as the Galerkin method is used. In this frame- 
work, the fields are expanded onto the Galerkin functions G. However the 
residuals R are expanded onto the standard basis. In the Galerkin method, 
the test functions are the Galerkin ones. One then gets as many equations 
as the number of Galerkin functions, formally written as {R, G) = 0. One 
can show that such equations are linear combinations of the coefficients of 
R which depend on the Galerkin basis used. As for the r-method, some 
equations can be relaxed to enforce additional conditions. 

To summarize, one can say that KADATH solves equations in the coeffi- 
cient space. Standard r-method is used to impose appropriate matching and 
boundary conditions. When regularity conditions must be enforced, this is 
done by means of a Galerkin method. 

4-2. Equations for the fields 

The discretization of the field equations by means of the r and Galerkin 
methods are implemented by the class Equation and its derived classes. 
Given the knowledge of the geometry and the type of equation, those objects 
are able to produce the appropriate number of discrete residual equations. 

The most widely used derived class is called Eq_inside. It deals with 
equations that must be solved inside a given domain, for instance a Pois- 
son equation of the type A/ — S = 0. For each type of domain, the r or 
Galerkin residual equations can be generated. Depending on the number of 
boundaries, several equations are relaxed in order to impose matching and 
boundary conditions. By default Eq_inside assumes that the equation to 
be solved contains second order derivative of the fields. For first order equa- 
tions (resp. zeroth order) one would use the related class Eq_one_side (resp. 
Eq_full). 

Boundary conditions are encoded in the class Eq_bc which must be sup- 
plied with the domain and the type of boundary, r or Galerkin method 



20 



of appropriate dimensions are constructed by this class. For instance, the 
boundary conditions on a sphere are written with a two-dimensional Galerkin 
method with respect to the angles {9, ip), in order to account for the regularity 
conditions on the axis. 

The class Eqjnatching is used to impose the matching of quantities across 
the boundary between two domains. This is done in the coefficient space and 
so must be used only when the basis relative to the surface is the same on 
both side of the boundary. This may seem restrictive but it covers most of 
the geometries implemented in KADATH , at least when the same resolution is 
used in each domain. This is the case for all the boundaries of the spherical 
geometry but not for the boundary between the bispherical domains and the 
spherical compactified one (see Sec. 12.31) . In most of the cases the quantity 
to be matched is simply an unknown field or its normal derivative. Never- 
theless it is possible to impose the matching of different quantities across the 
boundary. This would be useful in the case where different variables are used 
in different domains (see an example in Sec. \7.2\i . 

When the basis are not the same across the boundaries, because of differ- 
ent geometries or different resolutions, one must use the class Eq_matching_non_std 
that performs the matching in the configuration space. In this case, given the 
fact that the collocation points are usually different on the two sides of the 
boundary, one must choose the set of points at which the matching conditions 
are written. It is up to the user to make sure that this choice is consistent 
and leads to a global problem that is not over or under-determined. 

4-3. Integral equations 

KADATH also enables to impose global conditions on the fields. One could 
prescribe the value of the total energy content of space or the integral of 
a given field on some surface, for instance. If the associated equations are 
computed in terms of the coefficients of the fields, they are not functional 
equations and do not require the machinery of the weighted residual method. 
They simply translate into additional conditions that must be fulfilled by the 
fields. 

In order to ensure that the full system contains the same number of 
unknowns than conditions, integral equations are generally associated with 
global unknowns. This is for instance the case of black holes, where the local 
rotation rate is an unknown that is constrained by demanding that a global 
quantity, the spin, takes a given value (see Sec. 17. 3p . Another example is 
found in the critic solution case of Sec. 17. 2^ where the period of the solution 
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is sought by imposing that one of the modes of one of the fields vanishes at 
the center. 



5. Solving the system 

5.1. Newton-Raphson iteration 

It follows from the techniques described in Sec. [3] and HI that the solver 
translates a system of functional equations into a finite set of algebraic ones. 
In general, those equations are non-linear and KADATH relies on the well- 
known Newton-Raphson technique to find a solution (as in jo], [13] )• Let us 
consider a set of unknowns, formally denoted by the vector u. The system 
of equation can be written as / (m) = 0. Obviously the vectors u and / must 
have the same size. 

The solution is sought by iteration, starting from an initial guess iP. Let 
us denote the approximation of the solution found after n**^ iterations. 
The Newton-Raphson scheme proceeds as follows. First one computes the 
vector /" = / (m"). If /" is small in some appropriate sense (for instance the 
maximum norm is smaller than a given threshold), then is good enough 
solution. If not, one needs to compute the Jacobian matrix J" of the system 
defined as: 

dfi_ 

duj 



= ^ m (19) 



where fi and uj denote to the ith and jth components of / and u respectively. 
One row of the Jacobian corresponds to one equation, and one column to the 
derivation with respect to one of the unknowns. The matrix is computed at 
the position of the current solution it^ and so must be recalculated at each 
iteration. 

The linear system 

J"X" = r (20) 

must then be solved and the next approximation of the solution is given by 
{[ri+i _ ^ _ j^n rpj_^g method is known to converge rapidly to the solution, 
at least if the initial guess is good enough. For linear systems, the Newton- 
Raphson method finds the solution in one iteration step. 

5.2. Automatic differentiation 

For simple problems it can be possible to explicitly derive an analytic 
expression of the Jacobian matrix from the knowledge of the equations /. 
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This can be done by making use of some software that can do symbohc 
computations. Such analytic expression can then be fed to the code and 
used to compute the Jacobian. However, in practice, that can be tricky 
when the equations are complicated, especially if their explicit expressions 
in terms of the unknowns are intricate (this is especially true for Sec. 17. 4[ 
where an explicit expression of the equations in terms of the spatial metric 
would be almost impracticable). 

In the context of KADATH , the computation of the Jacobian is done nu- 
merically by a technique of automatic differentiation based on the notion of 
dual numbers. Each quantity x is supplemented by its infinitesimal variation 
6x. The new object is usually denoted as {x,6x). The arithmetic of those 
dual objects is then implemented in order to accommodate the usual rules 
for differentiation. For instance one gets: 

{x,5x) + {y,5y) = {x + y,5x + 6y) (21) 
{x,6x) X {y,5y) = {xy,x6y + y6x) (22) 

= (v^'^)- (23) 

Let us consider a set of values u of our unknowns and let us supplemented 
it with a set of infinitesimal variations 6u. Using the extended algebra one 
can apply the system of equations / to {u, 6u) . This would result in the 

following extended object (^f (u) ,6f {u,6u)^. The first part f (u) is the 

usual application of the system of equations to u and would be used to get 
the term /" in the Newton- Raphson iteration (see Sec. 15. ip . 

One can show that the second part Sf is a vector that is the product of 
the Jacobian with the vector 6u so that: 

6f{u,6u) = {J{u))x6u. (24) 

The use of the dual numbers enables to automatically compute the product 
of the Jacobian with any vector. In order to compute the whole Jacobian 
matrix, one then proceeds as follows. First the unknowns are supplemented 
with a variation vector that is zero except at the i^^ component which is 
set to one. The system of equation is applied, in its extended form, to this 
object. The variation of the result is then, by construction, the i^^ column 
of the Jacobian. By taking all the possible values of i, the whole matrix can 
be generated. 
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5.3. Parallelization and linear system solvers 

The size of the hnear system f l20|) is the same as the total number of un- 
knowns. If one considers Nf scalar fields, in domains of d dimensions and 
if N is the number of coefficients of the spectral expansion in each dimension, 
then the Jacobian is an (m x m) matrix with m ^ NfN^N'^. This number 
can be quite large, especially in three dimensions. For instance, for Nf = 5, 
Nd = Q, d = 3 and = 21, which are big but not huge numbers, one would 
get m > 250, 000. The resulting matrix would represent more than 500 GB 
of data. Such amount of data would be difficult to store on a single processor, 
without mentioning the number of operations required for the inversion. 

There are several ways to overcome the difficulty of dealing with such big 
matrices. One solution is to use an iterative technique to solve the system 
f l20|) (see [i^] for a general introduction to the iterative techniques and ji, 11 
for some applications). Doing so, the solution of the linear system is sought 
in a loop that does not require the storage of the full matrix J. It is only 
needed to be able to compute the product Jx for any vector x. If this is 
exactly the information that is available from our automatic differentiation 
technique (see Sec. 15. 2p . KADATH , in its current state, does not use this 
kind of algorithms. Indeed, the tests conducted with the iterative solvers, in 
the context of KADATH , showed a lack of stability, meaning that the conver- 
gence was achieved only for some problems and only for moderate degrees 
of freedom. The reason for that must probably be sought in the precondi- 
tioning step. Indeed, the success of the iterative algorithms relies on the 
knowledge of an efficient preconditioning matrix that approximates J~^. If 
some preconditioning techniques have been successfully used in the context 
of spectral methods like in ji, 11|, they are only applicable when working in 
the configuration space, which is not the case of KADATH . Nevertheless, for 
very large problems, it is likely that the use of iterative techniques would 
be highly desirable. This is why their implementations will be an important 
axis of future development of the library. 

That being said, in its current state, KADATH relies on a direct method to 
invert the linear system. In order to deal with very large matrices, distributed 
algorithms are available. Typically, each processor has the knowledge of 
only some parts of the matrix. In this context, a parallelization of KADATH 
is straightforward and is done via MPI. Indeed, the Jacobian is computed 
columns by columns, each computation being independent of the other. Each 
processor computes and stores a manageable set of columns. The size of the 
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Jacobian is only limited by the total amount of memory available. This part 
of the computation is perfectly linear in terms of the number of processors. 
After having tested several libraries, it appeared that a parallel version of 



LAPACK called Scalapack |21|, was the best suited. However, before calling 
the function that computes the solution the linear system, a redistribution 
of the matrix amongst the various processors must be done. The reason for 
this is purely computational and its goal is to ensure a better behavior of 
Scalapack algorithms. Let us precise that Scalapack library does provide 
the function needed to do this redistribution. The LU decomposition of J 
itself does not require more memory than what is needed to distribute the 
Jacobian across the various processors. Concerning the computational time, 
it has been observed that the resolution of the linear system takes roughly 
as much time as what is needed to compute the Jacobian itself, at least for 
the cases exhibited in Sec. [3 

For small problems (typically 2-dimensional problems in low resolution), 
a sequential version of KADATH can be used, where the linear system is solved 
using the standard LAPACK library (22| . 



6. User interface 

If KADATH is a rather intricate tool, a great effort has been made in its 
design to render its use as simple as possible. This is especially true in the 
way the equations are passed to the solver. In this section the various steps 
in which KADATH should usually be called are briefly summarized. KADATH 
is written in C++ and it is probably best if the user has at least some 
knowledge of this language. What follows does not constitute an extensive 
documentation of the library but is intended to give a feeling of what a typical 
use of KADATH is. For more details, the would-be user is strongly advise to 
take a look at the codes that have be made available in KADATH repository 
(this is the case of the four examples discussed in Sec. [7]). 

The first step is to specify the geometry of space. This is done by calling 
the constructor of one of the derived classes of Space (see Sec. |2]). Various 
parameters describing the geometry are required at this point, like the num- 
ber of shells for a spherical space, or the various radii involved in a bispherical 
one. The constructor must also be supplied with the required resolution (i.e. 
the number of points) and the type of polynomials to be used (Legendre or 
Chebyshev) . 



25 



On the desired geometry, the user needs to define the fields of interest. 
They can be of various type, scalars, vectors, higher order tensors, metric 
tensors, etc. As the solution is sought by iteration, those objects are usually 
set to some initial guess. The choice of the initial values depends very much 
on the problem. It can come from reading data from another code or by 
finding appropriate analytic expressions. A crucial point of this step is to 
affect each field with its correct spectral basis of decomposition. For most 
cases, KADATH provides functions that do this automatically but the user 
should make sure that those functions are appropriate for each given problem. 
Failing to provide the right basis will either cause the solver to abort (if it 
needs to sum two quantities with incompatible basis for instance), or cause 
the appearance of Gibbs-like phenomenom (if a general function is expanded 
onto symmetric basis for instance). 

All the informations needed to solve a system are contained in the class 
System_of _eqs. It is constructed from the space of interest and can be sup- 
plied with the domains on which the system is to be solved, in case the 
values of fields are not defined on the whole space. The list of variables and 
constants are then passed to the system. By variables one means quantities 
that are unknowns and that must be solved for, whereas constants are quan- 
tities that can appear in the equations but have fixed values. Both type of 
objects can be either fields or numbers. Each variable and constant must be 
a quantity that has been initialized beforehand. Each of them is supplied 
with a character string that is the name by which it will by recognized in the 
equations. A special field is the metric one, when defined, for which addi- 
tional functionalities are available (tensorial indices manipulation, covariant 
derivative, etc.) 

The equations are passed to the System_of _eqs as character strings. 
Those strings are read by KADATH to generate the appropriate computa- 
tion rules. The formalism used is inspired by LateX p^], especially in 
the way indices are handled. Einstein's convention of summation on re- 
peated indices is used (for instance, in the 3-dimensional Cartesian case, 
ft* 9^ = fxQ^ + fyd^ + fzQ^)- Various reserved words are available to encode 
some functions like sqrt that stands for the square root, dn for the normal 
derivative with respect to one surface, or R for the Ricci scalar of a given met- 
ric. When a quantity appears many times in the equations, it can be made 
into a definition. Apart from the fact that it can simplify the writing of the 
equations, it has the advantage that definitions are computed only when the 
value of the variables change, and not each time they are encountered, thus 
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resulting in a faster computation. When setting the equations, the user must 
not only make sure they are correct but also that they give raise to as many 
conditions as the number of unknowns, in order to get an invertible system. 
This can be tricky at times, especially when different resolutions are used 
in different regions of space. In some standard cases, member functions of 
the Space class provide means to implement equations in the whole space 
instead of passing them domain by domain. 

Finally a solution of the system is sought, calling the member function 
do_newton, the required accuracy being passed as a parameter. At this point 
KADATH verifies that the number of unknowns is consistent with the number 
of equations before entering the Newton-Raphson loop. Once convergence is 
achieved, the various data can be saved, manipulated, results printed, etc. 

7. Test problems 

In this section, four different problems are presented. The main goal 
is to illustrate the ability of KADATH to deal with different and complicated 
situations. The examples have been chosen to show various and different as- 
pects of KADATH . Nevertheless, they are far from being merely test problems. 
They all related to very contemporary physics. That being said, this paper 
is not the place to extensively discuss physics and so explanations about this 
are kept to the minimal level required for comprehension. Given the cur- 
rent limitations of KADATH , the four problems have in common the fact that 
they all are boundary value problems on a fixed and known geometry. The 
results shown below were obtained by using Chebyshev polynomials. The re- 
sults were roughly the same with Legendre polynomials, even if the observed 
convergence seems to be slightly slower in this case. 

Computations were conducted on a cluster of 10 quadcore bi-opterons 
at 2.5 GHz. Each node has 16 Go of RAM and the cluster runs under 
Linux. Under those conditions, the biggest Jacobian invertible has a size just 
below 100, 000 and it takes 4 h to do the inversion. Most of the codes used 
to compute the following results are available from the Codes/Par_version 
repository of KADATH . 

7.1. Vortons 

In classical field theory, some closed vortex loops can exist. Such loops 
are stabilized either by twisting them, in which case one talks about knots, 
or by centrifugal force when the vortex is spinning. This last case is called a 
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vorton (see 23| for a review on this type of objects). In this section, vortons 



are computed in the so-called Wittens theory [2J]. In this model, super- 
conducting cosmic strings are known to exist. They are, basically, straight 
configuration of the fields. The idea behind vortons is to try to cut a piece of 
a string, make it into a loop and try to stabilize it by giving it a spin. If this 



idea dates back to |25l. |26||. it is only recently that the such solutions were 



explicitly computed in [23j. 

From a mathematical point of view, vortons are described by two complex 
scalar fields and o. The vorton is computed by demanding that a takes 
the particular form 

a = Z exp (my9 + cut)] (25) 

where Z is a real function, m is the azimuthal winding number that con- 
straints the topology of the vorton, ip is the angle around the axis of the loop 
and u is the angular velocity. No particular ansatz is assumed for the other 
field that is described by its real and imaginary parts: = X + iY. 

The geometry of the solution is such that the fields are axisymmetric. 
Moreover, X and Z are symmetric with respect to the plane z = and Y 
is antisymmetric. Let us note that from the point of view of KADATH , Z is 
not a scalar field but rather the harmonic m of a scalar field (i.e. Z cos {rmp) 
is the real scalar field). This has to be taken into account when setting the 
appropriate spectral basis. This also changes the regularity conditions on the 
axis, where Z must vanish, which would not be the case for a real scalar field 
(see Sec. 13.41 for more details). Given the geometry, the problem is defined 
on a polar space by the class Space_polar (see Sec. 12. 2p . 

Under those assumptions, the three unknown fields obey a set of three 
elliptic equations: 



AX = (^^{X^ + Y'-1)+^Z^^X 
AY = (^^(X^ + Y'-1)+^Z^^Y 
A^Z = (^^{Z'-vl)+l{X' + Y^)-co'^Z 



(26) 
(27) 
(28) 



where A^ is the Laplacian of the mth-component of a scalar field, so that 

AmZ = AZ 2~- This expression shows that Z must vanish on the 

sin 6 

z-axis. At spatial infinity the fields are such that X = 1, Y = and Z = 0. 
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The equations are passed as such to KADATH . Being second order dif- 
ferential equations, matching of the fields and of their normal derivatives is 
performed at the boundaries between the various domains. The main diffi- 
culty in finding solutions lies in the fact that the system involves many free 
parameters X^, X^, 7, ?7o-, 1^ and m. It is believed that solutions only exist 
in some parts of this huge parameter space. Another difficulty is related to 
the existence of trivial solutions (like X = 1, Y = and Z = in the whole 
space). If one does not start the solver from a good enough initial guess, 
then the interesting solutions will be missed. For those reasons, it is quite 
an achievement that the authors of 23| did identify an appropriate region of 
the parameter space and were able to compute the associated vortons. 

The results presented in this section aim at reproducing some of the 



results obtained in [23[ in the case m = 2. The other parameters cirG clS 
follows A<^ = 41.12, Xa = 40, r/o- = 1 and 7 = 22.3. One configuration 
computed in [23| is used as an initial guess and a sequence of configurations 
rotating at different speeds is then constructed by slowly varying u. The 
precision reached by the code is asserted by computing the deviation of the 
solutions with respect to some virial theorem. The virial error is defined as 

1 - 3 ~ -^0) (ggg (6.159) of [23]). A^, Eq and E2 are integrals of 

-C/9 



the fields which explicit expressions are given by Eqs. (6.154) of 23|. The 



virial error is shown in the left panel of Fig. [21 as a function of u, for three 
different resolutions. The convergence of the virial error is rapid, as expected 
with spectral method. The error is below 10~^ in the whole range of u, for 
the highest resolution. 

The right panel of Fig. |2] shows the total energy E an the Noether charge 
Q of the vortons, as a function of u, as defined by Eqs. (6.155) and (6.157) of 



23|. The results are consistent with Fig. 30 of |23|]. For the low values of cu, 
one is limited by the fact that the vortons are bigger and bigger. One would 
need to use more domains to get to smaller u. At the higher end, there seem 
to be a transition from the vorton solutions to solutions for which Y = 
in the whole space. This transition prevents the code from reaching as high 



values of 00 as in j23|], for a yet unknown reason. Figure [3] shows contours 
of constant values for the fields for u = 0.85. Let us finally mention that 
vortons with m = 1 were also successfully computed. From the numerical 
point of view, they mainly differ from the m = 2 case from the fact that the 
basis of decomposition for Z is different. 
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Figure 2: Sequence of vortons with different uj. The other parameters are m = 2, 
= 41.12, Act = 40, rja = 1 and 7 = 22.3. The left panel shows the virial error for three 
different resolutions. The curves are labeled by the number of points in each dimension. 
The right panel shows the total energy E and Nocther charge Q along the same sequence. 



7.2. Critical collapse 

Critical collapse was first observed in 27| when computing the evolution of 
a spherically symmetric massless scalar field in general relativity. In the weak 
field regime the field disperses at infinity whereas, in the strong field regime, a 
black hole is formed via gravitational collapse. Just at the threshold between 
those two cases, another type of solution appears. It is a naked singularity 
and is called the critical solution. It has many interesting properties and we 
refer the reader to [28] for a review on this subject. The structure of the 
critical solution for a scalar field collapse has been extensively studied, for 
instance in 29|, |30|. This section is devoted to the computation of the critical 
solution by KADATH . 

As stated in |3(|], the critical spacetime can be found as the solution of a 
2-dimensional non-linear problem. Coordinates can be chosen such that one 
spatial variable x goes from to 1 and one periodic variable r goes from to 
27r. r is the time coordinate. From the geometrical point of view, the space 
of interest is a cylinder. The relevant class in KADATH is called Space_critic 
and is described in detail in Sec. 12.41 In order to take into account some 
symmetries at x = 0, this space is separated into two domains. 

The solution involves four fields, two describing the matter U and V and, 
two for the gravitational field, / and a. The fields have some symmetries 
with respect to r = vr so that one can restrict r to [0,7r[. More precisely 
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Figure 3: Fields X, Y and Z (from left to right) for uj — 0.85. The other parameters are 
the same as in Fig. [2] The continuous (resp. dashed) lines show constant positive (resp. 
negative) value of the fields. 



a (r + vr) = a (r) (idem for /) and f/ (r + vr) = —U (r) (idem for V). Some 
symmetries are also present at x = where the metric fields / and a are even 
function. Some combinations of the matter field are also even and finite at 
x = 0: U={V + U)/ (2x) and ^ = (V - [/) / {2x^). 

In order to take the symmetry at x = into account, different sets of 
unknowns are used in different regions of space. In the inner domain, one 
works with 11 and \l/ and in the outer one with U and V. Doing so, the various 
fields are expanded into the appropriate spectral basis of decomposition, as 
described in Sec. 13.2.31 

In the outer domain, the equations are given by 

= («'-!)/ (29) 
X (a-2) ^ = 1- {l + U^ + V^) (30) 

x{f + x)U, = f[{l-a')U + V]-'^xUr (31) 

x{f-x)V, = f[{l-a')V + U]+^xVr (32) 

where A is the period of the solution that will be discussed later. 
In the inner domain, one gets 

xU = (a'-l)/ (33) 
X (a-2) ^ = 1 - (1 + 2x^U^ + 2x^^^) (34) 

x{f-x')U^ = -U{f~x')+f{l-a'){fn + x'^) (35) 
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+/ (/n - x'^) + (/vi>_, + n 

x{f-x')^,, = -2^{f-x')+f{l-a'){f^ + U) (36) 
+f(Il-f^) + ^{fU^ + x'^,^). 

The equations are first order differential equations. The variable r being 
periodic, no boundary conditions is needed. The situation with respect to 
X is more complicated because some of the factors in front of the derivative 
vanish at the boundaries. The appropriate number of boundary conditions 
must be determined by a precise examination of the equations. 

Let us first turn to the variable a. On the side a; = 1 of the cylinder, 
Eq. fl30|) is not degenerated, so that it is treated in the standard way, that is 
standard first order r-method (see Sec. 14. ip . At the inner boundary however, 
Eq. flM|) is degenerated and becomes = 1. The equation has to be treated 
by a r-method of zeroth order (all the residual equations are kept). In a 
sense the equation is its own boundary condition. 

For the other metric field /, in the outer region, the situation is also the 
standard one. In the inner part, Eq. ( 15^ is also degenerated and becomes 
/ (a^ — 1) = 0. As previously seen, this condition is already accounted for by 
the equation for a at x = 0. This implies that the degeneracy is not effective 
and that Eq. f l33|) must be treated by a standard first order r-method. So 
the field / does require a boundary condition on one side of the cylinder. It 
is physically motivated and simply is / (x = 1) = 1. 

Given this value of / at the outer side, one of the equations for the matter 
fields equation ( l36l) is actually degenerated and equivalent to the regularity 



condition (20) of [30] (the sign in j30| is faulted by a typo). On the inner 



side equation ( 135|) is degenerated but equivalent to pU (1 — a^) = and is 
already accounted for by the equation for a. The other matter equation ( 136|) 
is also degenerated and gives a true condition /II— 3/^\l/+27r/A/n = 0. So, 
concerning the matter fields, no additional boundary conditions are needed 
for there are two degeneracy conditions (one on each side), for two fields. 

As one is dealing with first order equations, the fields themselves must 
be matched at the interface between the two domains. This study is slightly 
technical but it ensures that the constructed global system of equations is 
well-posed and invertible. The final situation is summarized in Tab. [2J Let 
us mention that such complicated setting (non-trivial matchings, different 
order in the r-method...) can be easily encoded in KADATH . 
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Table 2: Construction of a well-posed problem. The columns a; = and a; = 1 summarize 
the behavior of the equations at the boundaries. The columns labeled inner and outer 
give the order of the associated r-method. The matching conditions at the interface are 
also given. More details about this can be found in the core of the text. 



Field 


x = 


inner 


matching 


outer 


X = 1 


a 


dege. ^ = 


1 


Qth. 


= a+ 




non dege. 


f 


dege. 4^0^ = 


1 




r = /+ 


-|^st 


/ = 1 


matter 


n : dege. 


= 1 


For n : 


2xn- = v+ + 


f/+ 


U : 1"* 


U : non dege. 




^ : dege. 




^ : 0*^ 


2a;2^- = 1/+ - 


■U+ 


V : O'^* 


V : dege. 



Along with the field equations, there is also a global condition that con- 
straints the value of A and that demands that the mode corresponding to 
cos (2r) in the expansion of / vanishes at x = 0. This last condition falls 
into the category of the integral equations discussed in Sec. 14.31 It is be- 
lieved that the equations are very sensitive to the initial guess and that they 
will fail to converge to the critical one if one starts too far from it. This is 
the reason for using the solution computed by 

a 

starting point of the 

Newton-Raphson iteration. The errors that appear when passing data from 
one code to the other introduce enough numerical noise to make the loop do 
a few iterations. 

The precision of the code is assessed by comparing the value of the period 
A to the one given by Eq. (58) of [30]. The relative difference between the 
two is shown in the left panel of Fig. HJ as a function of the number of points 
in each dimension. The convergence is satisfactory even though it is difficult 
to state that the evanescent regime has been reached. Let us note that the 
convergence of the solution is significantly slower than for other problems, an 
accuracy of 10"*^ being reached with as many as 96 points in both dimensions. 
This is expected, given the strong gradients appearing in the fields. Those 
gradients are responsible for a rather low decay of the coefficients of the 



spectral expansions. This was already observed in [30| and can be seen in 
their Fig. 8 and 9. As an illustration, some fields at both sides of the cylinder 
are shown in the right panel of Fig. |H This plot is the equivalent of Fig. 4 



of 30|, except for the definition of the coordinate r. In [30\, the temporal 
coordinate r goes from to A, whereas in this work, r/A is used so that our 
coordinate goes from to vr, once the symmetry at half-period is taken into 



account. 
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Figure 4: The left panel shows the relative error in computing the period A. The exact 
value is assumed to be given by Eq. (58) of (s^l- The convergence is shown with respect 
to the number of points in each dimension. The right panel shows the values oi f {x — 0) 
(continuous line), ^ {x — 0) /300 (dotted line) and U {x — 1) (dash-dot line), as a function 
of T. 



7.3. Binary black holes 

There is a long history of works that aim at computing the structure 
of spacetimes that contain a system of binary black holes. Amongst the 
many reasons to study such systems, one will mention the fact that they 
are known to be good emitters of gravitational waves and so are one of the 



main target for the gravitational wave detectors currently in operation [31 



Binary systems have also a great interest in the context of galaxy formation. 
Indeed, it is believed that nowadays galaxies have been formed by successful 
mergers of smaller galaxies. When the merger occurs, the black holes present 
at the center of the smaller galaxies will become bounded and form a binary 



system [32 



Under the influence of gravitational radiation, two black holes will not 
remain on closed orbits but rather spiral towards each other, until they merge 
into a single object. This process occurs in the strong field regime of gravita- 
tion and must be described in the context of general relativity. Most of the 
simulations are performed in the 3-1-1 formalism where a splitting of space 



and time is introduced [331]. The main effect of this splitting is to separate 
Einstein's equations into two sets: (i) the constraint equations that do not 
involve time (ii) and the evolution equations that contains Dalembert type 
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operators. Doing so, the simulation of a binary system proceeds in two sep- 
arate steps. First, one needs to produce a initial configuration that must 
satisfy the constraint equations and that represents the physics of interest 
as accurately as possible. This is known as the initial data problem js^. In 
a second step the behavior of the fields at latter times is obtained by using 
the evolution equations. Let us mention that if the constraint equations are 
fulfilled at the initial time, they are verified at all time steps, if the evolution 
equations are used properly. 

In this section, one is interested in preparing an initial configuration that 
represents two black holes in quasi-circular orbit along the same lines as in 



35[. Even if true circular orbits cannot exist due to gravitational radiation, 
it is believed to be a rather good approximation, at least when the black holes 
are relatively far apart. Technically this is done by imposing the existence 



of an helical Killing vector (see for instance [36|, |37[ and references therein). 
Another approximation usually done in this context is the so-called conformal 
flatness approximation (see below for a precise definition). This has more to 
do with the substantial mathematical simplification that results than with 
any physical motivation. It was once believed that the conformal fiatness 
approximation would minimize the amount of spurious radiation in the data 
sets but it turned out not to be the case. Nevertheless, the approximation 
leads to a consistent mathematical problem, especially as far as the behaviors 



of the fields at infinity are concerned [37 



In the 3-1-1 formalism, the 4-dimensional metric g^^, that describes the 
system is expressed in terms of purely spatial quantities. More precisely the 
line element is given by 

ds^ = g^^dx^'dx" = - (A^2 _ jyijY;. j ^ 2Nidtdx' + 'yijdx'dx^ (37) 

where the Greek indices are 4-dimensional and run from to 3 and the 
Latin ones are purely spatial and go from 1 to 3. The 3+1 quantities that 
describe the geometry are one scalar the lapse, one vector A^' the shift 
(three independent functions) and one spatial metric 'jij (six independent 
functions). The conformal fiatness approximation states that the spatial 
metric relates to the fiat metric by a single scalar conformal factor More 
precisely one demands that 

1^J = (38) 

where fij is the fiat metric. This condition is only an approximation and it 
is not believed to be exactly true. 
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In this setting, the unknowns are the fields A^, \1/ and N'^. The associated 
equations are obtained from Einstein equations by enforcing stationarity (i.e. 
the time derivative are set to zero). It gives a set of five elliptic equations: 



DiD'N 



(39) 
(40) 
(41) 



where D is the covariant derivative with respect to the fiat metric, fij is 
also used to raise and lower indices of tensors. A^^ is the conformal extrinsic 
curvature tensor and represents the way the three metric is embedded in the 
4-dimensional geometry. From the mathematical point of view it is defined 
as being 



, - - . , - . - . . , . (42) 
2A^ V 3 - . ; \ ) 

The elliptic equations must be supplemented by appropriate boundary 
conditions. At infinity, it is demanded that fiat spacetime is recovered which 
implies that 



N 



= and \E' = 1 when r — )■ oo. 



(43) 



The inner boundary conditions are enforced on two spheres that represent 
the black holes themselves. They are basically obtained by demanding that 
those two spheres are apparent horizons in equilibrium (see 38| for a review) . 
Those inner boundary conditions enforce the presence and the physical prop- 
erties of the holes. On the sphere Sa {a = 1 or 2), one gets 



N\ 

a7 ^ 2^ 



Sa 

Sa 
Sa 



no 



^3 . . 



(44) 
(45) 

(46) 



Equation ( I44j) is just a choice of time coordinate on the spheres and uq is 
a constant (in this section one takes uq = 0.1). Equation fH5|) translates the 
fact that the spheres are apparent horizons, is the unit vector normal to 
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the sphere which reads [xajra-, Haifa, ^al'^a) in Cartesian coordinates centered 
on the sphere a. Equation f l46l) states that the spheres are in equihbrium and 
also contains information about the state of rotation of the holes, is the 
orbital velocity of the system and M* the constant vector (0,Xa,0) where 
Xa is the coordinate distance between the center of mass and the center of 
the hole a. This part of the boundary condition accounts for the orbital 
motion of the holes. is the local rotation rate of the black hole a and it is 
associated to the local vector = (— ?/a,Xa,0). This last part accounts for 
the spins of the objects. 

f2 and the two fia are global unknowns that must be constrained by 
additional equations, as discussed in Sec. 14.31 The proper value of the 
orbital velocity is obtained by demanding that two integral quantities, the 
Komar mass and the ADM mass are equal. This is closely linked to a virial 



theorem, as discussed in [37j.The two masses are defined as surface integrals 
at infinity 

MKomar = ^ / N ^S' (47) 

Madm = ^ / (48) 

In this work, the local rotation rates are determined by demanding that the 
black holes are not rotating and so that their spin Sa vanish. The individual 
spins can be defined by integrals on the spheres themselves 

Sa = ^(f ^'A,m]^dS^. (49) 

The total angular momentum is defined by a similar integral, taken at infinity, 
and reads 

J= — (( AijTn'dS^ (50) 

where m* = {—Y, X, 0). X and Y are the Cartesian coordinates with respect 
to the center of mass. 

This problem has been solved in the context of KADATH by using the bi- 
spherical space implemented by the class Space_bispheric (see Sec. 12. 3p . 
The input parameters are the radii of the spheres (taken equal in this partic- 
ular case) and their separation. A Cartesian tensorial basis of decomposition 
is used. The unknown are the fields A^, \1/ and A^' = {N^ , , N^) and the 
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global quantities Q and fi^. Aij is defined as being a definition in terms of 
the unknowns (given by Eq. fl42p ). The various vectors appearing in the 
equations, like the rrii are passed to KADATH as constants. The equations are 
given by Eqs. fl39t |40| |4T]) and the boundary conditions enforced on both 
the spheres and at infinity. Matching of the fields and their normal deriva- 
tives at the boundaries between the domains is imposed. The three integral 
equations constraining the global quantities are also passed to the solver. 

Contrary to the other cases presented in this paper, the binary black hole 
configurations are really three-dimensional. This implies that the computa- 
tional task is somewhat harder. Configurations with four different resolutions 
have been computed with 9, 11, 13 and 15 points, respectively, in all dimen- 
sions and in all the domains. Let us mention that the size of the Jacobian for 
the higher resolution is just less than 100,000. Configurations for five differ- 
ent separations are computed. It appears that the values of the fields depend 
strongly from the value of fl which must be obtained with a good accuracy. 
Such accuracy can be measured by looking at the convergence of Q as the 
function of the number of points as in Fig. O The convergence is shown by 
taking the difference between the value of Q found for given resolution and 
the one for the best available resolution (15 points in each dimension in this 
case) . 

The physical quantities are usually given in their adimensional form. This 
can be done by making use of the area mass of the black holes: m = ^y~A/TQ^T. 
A is the area of the hole which is given by an integral on the horizon (i.e. 
the sphere) : 

A = [ ^MS. (51) 



The area mass of the system is simply M = mi -|- m2 and is constant along a 
sequence of varying separation (see for instance 3^). It follows that M is the 
total mass of the system when the black holes are infinitely separated. One 
can then define the binding energy of the system as Eh = Madm ~ M. The 
reduced mass is defined in the usual manner and in the case of equal mass 
black holes is given by fi = M/4. The adimensional binding energy Eh/jj, 
and the adimensional total angular momentum J/Mfi are shown on Fig. [6l 
as a function of the adimensional orbital velocity MQ. The values obtained 
by KADATH are compared to the data from [s^l • The agreement is very good, 
especially considering that the variations of J (resp. Eb) are small compared 
to the values of J (resp Madm) itself. 
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Figure 5: Convergence of 17 as a function of resolution (compared to the best available 
value obtained with 15 points in each dimensions). The convergence is somewhat faster 
for small separations. 

Finally, contours of constant values of some fields are shown in Fig. [71 in 
the orbital plane (the plane z = 0) for a separation of 12. The locations of 
the spheres are clearly visible. 

7.4- Kerr problem 

In this section, one aims at recovering the exact solution for a single 



stationary rotating black hole, along the lines of |4C|. It is called the Kerr 
solution and is known to be analytic, at least for some choices of coordinates 
j4l|. This is however not the case for the formalism that is used here. The 
starting point is the same as the one used in the binary system case in Sec. 
17.31 An orthonormal spherical tensorial basis of decomposition is used. Doing 
so and for this problem, all the fields will independent of thus recovering 
the fact that the solution is axisymmetric. 

In this section, and as opposed to what is done in Sec. I7.3[ one aims 
at recovering an exact solution and so one needs to remove the conformal 
flatness approximation. Indeed, even for a single black hole, it is known that 



39 



0.12 0.14 0.16 0.08 0.1 0.12 

Mii Ma 



Figure 6: Binding energy (left panel) and total angular momentum (right panel), as a 
function of Mfl. The circles denote the values obtained by KADATH and the solid curves 
are data taken from 1351. 



their exists no choice of coordinates for which the conformal metric is flat 



42]. That being said, the unknowns are the same as in Sec. 17.31 : two scalars 
the lapse N and the conformal factor and one vector, the shift N\ to 
which one must add the conformal spatial metric 7jj itself. The metric is a 
rank-2 symmetric tensor and so has six different components. By definition 
of the conformal factor, is such that its determinant is one. As will be 
seen in the following, the conformal metric still contains some gauge degrees 
of freedom. 

The constraint equations along with the trace of the evolution ones can 
be written as 

DiD'N = -2 ^^ +N^^A,,A'^ (52) 
= (53) 

b^^ 

D^A, = (54) 
where A^^ is the extrinsic curvature tensor which is given by 

A'^ = ^ (^b'N^ + b^N' - '^bkN'^f^^ (55) 
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Figure 7: Values of the lapse function (left panel), the conformal factor (center panel) 
and the component Ny of the shift (right panel), in the orbital plane z = 0. The coordi- 
nate separation is 12. The continuous (resp. dashed) lines show constant positive (resp. 
negative) value of the fields. 



where D denotes the covariant derivative associated to 7jj and R is its scalar 
Ricci tensor. Equations f l52|) to fl55l) are the equivalent of Eqs. fl39l) to fH2|) . 
in the non-conformally fiat case. 

At infinity one demands that fiat spacetime is recovered which gives N = 
1, \E' = 1 and iV* = 0, as in Sec. 17.31 The inner boundary conditions are 
obtained in the same manner as in the binary case except that they must be 
written with a general (non-fiat) spatial three metric. One then gets, on the 
sphere S 



AS' 



no 



\]/2 



(56) 
(57) 

(58) 



Equations fl56l) to fl58l) are the equivalent of Eqs. (jHj) to fH6l) . in the non- 
conformally fiat case and for a single blacic hole, s is the unit outgoing 
normal to the sphere, with respect to the conformal metric and Q the rotation 
parameter of the hole. 

The evolution equations, written in the stationary case, can be seen as 
equations for the spatial metric and they are 



DiDjN - 2- 



DiND.m 



- 2- 



+ 27, 



6iV^^ (59) 
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+2N^^ + 2iV%- + 2iV%-^^^ + 2iV^M,fc4 

where is the Ricci curvature of 7jj. It is the curvature that contains 
the second order operators in terms of the metric (i.e. the expression of Rij 
involves terms that look like A (7ij)). 

Equations ( 159|) are obviously symmetric and so represent six components. 
However, if used as such, they do not lead to a well-posed problem because 
they are not all independent. The construction of a well-behaved problem is 
rather tricky in this case and is related to both the inner conditions for the 
spatial metric and the choice of gauge. 

The problem of knowing what boundary conditions must be enforced on 
the metric on the horizon is a very contemporary one and several proposals 
have been made in the literature. However only one has been successfully 
applied so far and it is known as the no-boundary treatment [40)]. The basic 
idea is to assume that Eqs. flS^ are somewhat degenerate on the horizon 
(like some equations of Sec. I7.2p and so that they require less boundary 
conditions. If the degeneracy is not strictly demonstrated, there are some 
good reasons to beheve that it is true, at least at the first order (see section 
IV-B of 40]). However the code used in 4^ does not use Eqs. ^9\} as such 



but a variation of them that aims at imposing that the spatial metric fulfills 
a particular gauge known as the Dirac gauge. This gauge can be written as 

Vif^ = 0, (60) 

where V denotes the covariant derivative with respect to the fiat metric. The 
gauge is used to simplify the explicit expression of Rij. 

When the full set of Eqs. (159|) is used, a full no-boundary treatment 
leads to a non-invertible system of equations because the gauge has to be 
imposed at some point. The idea is to use Eq. (1^ . or at least some part of 
it, as a boundary condition for the metric. From the numerical experiments 
conducted with KADATH , it is found that a way of recovering the Kerr space- 
time is to use only the component 6 of Eq. f l60|) as a boundary condition 
for the component {6,r) of Eqs. fl59l) . The other components of Eqs. fl59l) 
are treated as degenerate without any boundary condition. It is somewhat 
surprising to observe that imposing only one component of the Dirac gauge 
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and only on the horizon is sufficient to obtain the full gauge in whole space, 
as will be checked a posteriori. Finally one has to note that the conformal 
metric 7jj must be, by definition of such that its determinant is unity: 
det (7) = 1. This is an additional equation on the components of the metric 
and it is used in place of the (v?, v^) component of Eqs. fl59l) . 

To summarize the situation, the equations that one uses when solving 
this problem in the context of KADATH far as the metric is concerned: 

• the (r, r), (r, (9,9) and {9,{p) components of Eq. fl59|) . without any 
boundary conditions near the horizon (i.e. a r-method of first order is 
used in the inner domain). 

• the (r, 9) component of Eq. fl5^ is used in the standard way, with the 
component 9 of Eq. ( l60i) as an inner boundary condition. 

• the condition on the determinant of the metric provides the last equa- 
tion. As it contains no derivative, it is solved by a r-method of zeroth 
order. 

At spatial infinity, one demands that the metric goes to the fiat one, which is 
possible only in the single black hole case (for a binary system, gravitational 
waves would probably forbid this). 

The reason why the setting described above leads to a system of equations 
that behaves correctly and why it is the only combination that does so, at 
least in the context of this paper, is currently not known. Further studies 
and deeper mathematical understanding of the system is probably required 
but beyond the scope of this work. 

Using a spherical space, those equations have been solved using four dif- 
ferent resolutions (9, 11, 13 and 17 points in each dimension). Sequences of 
black holes rotating at different values of fl have been obtained. One single 
spherical shell is used and the value of Uq is set to 0.5 in all the cases (see 
Eq. flS^ ). The results are presented as a function of the usual Kerr parame- 
ter J/M^ = a/M where J is the total angular momentum and M the ADM 
mass of the system. Both quantities are computed by surface integrals at 
infinity (see Sec. 17. 3p . The value of the Kerr parameter as a function of the 
adimensional quantity MQ is plotted in Fig. [SI A value of a/M ^ 0.91 is 
reached. This value is lower for lower resolutions as can be seen in Fig. [HI 

Two error diagnostics are shown in Fig. [91 On the left panel the residual 
error on the full set of equations is plotted. By full set, one means not only 
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Figure 8: Value of the Kerr parameter as a function Affi for the high resolution configu- 
rations. The values for lower resolutions lies almost exactly on the same curve. 



the equations used to solve the system but also the ones that were forgotten. 
In particular one can think of the (y?, (f) component of Eq. ( 159|) or the Dirac 
Gauge in whole space. This is a very strong test and it shows, as already 
stated, that the Dirac gauge is enforced in the whole space by only its 9 
component on the horizon. The right panel of Fig. |9] shows the relative 
difference between the ADM and Komar masses (as defined by Eqs. (H7I) 
and (H8ll ). Like in the binary case, those two quantities should be equal. 
However, it is never directly enforced in the single black hole case and so 
gives a check of the precision of the code. As stated in Sec J7.3[ this is related 
to a virial theorem 37|. Both plots in Fig. |9]are shown as a function oi a/M 
and for the four different resolutions. The errors decay as the resolution 
increases and are higher for higher values of the Kerr parameter. This is 
probably related to the fact that the fields have stronger gradients when 
a/M increases and so need more points to be accurately described. This 
effect is also observed in 40|. Some examples of the fields, for a/M ^ 0.91 
are shown in Fig. [101 
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Figure 9: Precision reached by the code as a function of a/M, for the four different 
resolutions. The left panel shows the maximal error on the full set of equations and the 
right one the relative difference between the Komar and ADM masses. The curves are 
labeled by the number of points in each dimensions. 

8. Perspectives 

This presentation of the KADATH hbrary was intended to convince the 
physicists that it can be a valuable tool in the study of a wide class of 
problems where partial differential equations are involved. Several domain 
decompositions have been presented (Sec. [2]) and the associated spectral 
basis exhibited (Sec. [3]). The discretization of field equations by means of 
the r and Galerkin methods have been discussed in Sec. H] and the solution 
of the resulting non-linear system explained in Sec. |5l The construction of a 
code that uses the KADATH library is briefly sketched in Sec. |6l 

The resolution of four different problems has also been presented. Solu- 
tion called vortons have been computed in Sec. 17. using a standard KADATH 
setting with spherical domains. In Sec. 17.21 the computation of a critic solu- 
tion in the context of core collapse has been presented. This case illustrates 
the use of specialized domains and different variables in different regions of 
space. Section 17.31 explained the study of spacetimes containing two black 
holes. It makes use of the bispherical coordinates implemented in KADATH . 
The ability of the library to deal with intricate sets of equations is clearly 
illustrated by Sec. 17.41 where the Kerr black hole is recovered, in the S-f-l 
formalism of general relativity. 
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Figure 10: Values of the lapse function (left panel), the component iV of the shift (center 
panel) and of the component 7'''" of the spatial metric. The associated value of a/M is 
about 0.91. The continuous (resp. dashed) lines show constant positive (resp. negative) 
value of the fields. 

Even if KADATH is currently in a production state, where it can be used 
to study new and interesting physics, there are still some developments that 
come to mind. One can for instance think about implementing new geome- 
tries like true cylindrical coordinates (as opposed to the speciahzed ones used 
in Sec. I7.2p or deformed spheres. Such spheres have been extensively used in 
the LORENE library to match the physical surfaces of deformed objects, like 
rotating neutron stars. This case is in fact a situation where the surface of 
the domains in not known beforehand but must be determined by the code. 
To do so, there will be need to have some additional unknowns that would 
describe the shape of the various domains. Those unknowns would be associ- 
ated to additional equations constraining the shapes. In the case of neutron 
stars, one would for instance demand that the surface of the object is the 
surface on which the specific enthalpy vanishes. Given those techniques have 
already been successfully applied in many codes, they should work properly 
with KADATH . Geometries with more than three spatial dimensions could 
also be studied, in order to simulate things like brane worlds. Needless to 
say that in this last case, the resulting size of the system would be quite big 
and possibly difficult to handle. 

It would be also interesting to test the behavior of KADATH on much bigger 
clusters than the one used so far. Several such clusters (with much more 
than 1000 nodes) are available to the scientific community and there is hope 
to install KADATH on such machines. Another way to deal with very big 
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matrices would be to use iterative techniques to invert the Jacobian. As 
stated in Sec. 15. 3[ this solution is currently not retained in KADATH because of 
several difficulties (lack of generality and absence of simple preconditionning 
techniques). Nevertheless, a much more detailed study of those algorithms 
should be undertaken. A successful use of the iterative techniques would 
potentially result in a dramatic decrease of the resources needed to handle 
the Jacobian, thus making the library easier to use, especially in the three- 
dimensional case. 

A major extension of KADATH would be the possibility to perform time 
evolutions. Even when spectral methods are used to discretize space, time is 
usually dealt with by means of finite difference schemes. One of the reason 
for this lies in the fact that the interval of interesting time is usually not 
known in advance. Also the discretization with respect to time is usually 
believed to lead to smaller errors that the spatial ones, thus making the use 
of temporal spectral methods less appealing. If those mixed schemes (finite 
difference in time and spectral in space), could be almost directly used with 
KADATH , it is not clear what this will improve with respect to already existing 
codes. Indeed, for free evolutions (like in jisi]) for which the only equations 
solved are hyperbolic ones, the whole KADATH machinery would not be used 
but to compute the sources. On the other hand, for constrained evolutions 



(like in |4J|), one would need to solve elliptic equations at each time-step, 
which would probably be impracticable for three-dimensional problems with 
KADATH . The library however, could be very useful in investigating schemes 
that are fully spectral, in time as in space. Results of this type are very 



few and are only available for very simple systems (see 45| for an example 
with only one spatial dimension). For such fully spectral codes, time is 
just another coordinate that is treated in the same manner as the spatial 
dimensions. If this could in theory be solved by KADATH , there is a long road 
to compute complicated system evolutions in this framework (like the binary 
black holes for instance). Preliminary tests, using systems with symmetries 
(like rotating single objects or pure gravitational waves), would be in order 
to check whether such computations are doable or not. 

Obviously all the work done would be quite pointless without applications 
to real physics problems and it is the main axis of future work. There are 
plans to apply the library to the computation of solutions in gauge field 
theories, like monopoles with or without gravity, vortons in more complicated 
theories than the one exposed in Sec. 17.11 In this area of physics, the field 
of application is huge. KADATH has been designed especially with this kind 
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of things in mind and it is believed that it will be very efficient in findings 
those solutions. 

Another area of application is the one of general relativity, especially in 
the context of compact objects, black holes or neutrons stars. One of the 
main application would be the computation of binary black holes configu- 
rations without the conformal fiatness approximation. In a sense it would 
be a merging of Sec. 17.31 and 17.41 But there are some conceptual difficul- 
ties coming from the presence of outgoing gravitational waves in this case. 
The applications of KADATH to cases containing neutron stars would also be 
interesting. If the gravitational field is usually weaker than for black holes 
and so easier to handle, the inclusion of matter can cause many additional 
difficulties. 

In the hopefully long life of KADATH , there is hope that it will also be 
applied to problems that are not yet known to the author. This would fully 
test the modularity of KADATH . Should this be successful, then KADATH would 
have reached its main objective. 
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